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ABSTRACT 

We investigate the properties of the Lyman alpha forest as predicted by numerical 
simulations for a range of currently viable cosmological models. This is done in order 



j^ I to understand the dependencies of the forest on cosmological parameters. Focusing on 

the redshift range from two to four, we show that: (1) most of the evolution in the 
Lj distributions of optical depth, flux and column density can be understood by simple 

;_j ■ scaling relations, (2) the shape of optical depth distribution is a sensitive probe of the 



amplitude of density fluctuations on scales of a few hundred kpc, (3) the mean of the 
b distribution (a measure of the width of the absorption lines) is also very sensitive to 
fluctuations on these scales, and decreases as they increase. We perform a preliminary 
comparison to observations, where available. A number of other properties are also 
examined, including the evolution in the number of lines, the two-point flux distribution 
and the Hell opacity. 

Subject headings: cosmology: large-scale structure, intergalactic medium, methods: nu- 
merical, quasars: absorption lines 
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1. Introduction 

Several numerical simulations of the Lya forest in cold dark matter (CDM) dominated cos- 
mologies have been performed in recent years and compared with observations (Cen et al. 1994; 
Zhang, Anninos, & Norman 1995; Hernquist et al. 1996; Miralda-Escude et al. 1996; Zhang 
et al. 1996; Dave et al. 1997; Bond & Wadsley 1997; Zhang et al. 1998). Remarkably ah the sim- 
ulations have been able to reproduce the measured neutral hydrogen column density distribution, 
the size of the absorbers (Charlton et al. 1997), and the line number evolution reasonably well, 
despite the differences in the cosmological models used: Cen et al. adopt a ACDM model; Zhang 
et al. investigate sCDM models with both an unbiased and a cluster scale normalization; Hernquist 
et al. (1996) evolve an sCDM model with a cluster scale normalization. The distribution of Doppler 
parameters has fared somewhat less well: the predicted distribution peaks toward lower values than 
observed when the simulations are performed with adequate resolution (Bryan et al. 1998; The- 
uns et al. 1998). Nonetheless, the generally good agreement with observations of the Lya forest 
suggests that the models are capturing the essential physical properties of the absorbers. This has 
prompted recent work by Croft et al. (1998) aimed at using flux statistics of the observational 
data to extract the fluctuation spectrum of the underlying cosmology. We are thus encouraged to 
investigate the possibility that differences in the statistical properties of the Lya forest predicted 
by different cosmological models may provide a means of testing the models. 

The objective of this paper is to compare the Lya forest statistics derived from simulations 
in different cosmological models and to investigate what key properties of the cosmological models 
control a given statistic. We present results from nine numerical simulations using five different 
background cosmological models, three of which are flat with no cosmological constant, one is open, 
and one is flat with a nonzero cosmological constant. For five of the simulations, which we will refer 
to as the model comparison study, the parameters of the cosmological models have been selected 
by their ability to match the local or low redshift observations, although all of these models except 
the standard cold dark matter (sCDM) model are also consistent with COBE measurements of the 
cosmic microwave background. A tilted CDM model is further designed to match COBE constraints 
on the normalization of the power spectrum on large scales. In the remaining four simulations we 
keep the underlying cosmology fixed (sCDM) while varying the normalization of the fluctuation 
power spectrum in order to clarify the dependence of the Lya statistics on this parameter. The 
radiation field is normalized to the absorption properties of the Lya forest as measured at high 
redshifts. While the emphasis in this paper is on comparing cosmological models, we also test how 
well the models are doing using the tabulated statistics of the Lya forest as determined primarily 
by Kim et al. (1997) for several QSO lines-of-sight. A more complete comparison with existing 
data from several observational groups will be presented in Meiksin et al. (1999). 

The paper is organized in the following way. In §^ we describe the cosmological models and 
simulation technique. In §y we investigate the model differences, power dependence, and redshift 
evolution in the raw opacity data as characterized by nonparametric statistics of the flux and 
optical depth. In §Q we present a line analysis of the spectra generated by the various simulations 



focusing on the column density distribution and line number evolution statistics. In §^ we discuss 
the Doppler b parameter distributions and related nonparametric statistics, and in §^ we present 
model predictions for He ii absorption. We summarize our results in §^. 



2. The Models and Simulations 

All the model background spacetimes we consider are in the context of Cold Dark Matter 
(CDM) dominated cosmologies. We examine the following five models: a standard critical density 
flat CDM model (sCDM), a flat CDM model with a nonvanishing cosmological constant (ACDM ), 
a topologically open CDM model (OCDM), the standard CDM model but with the power spectrum 
of the density perturbations tilted (tCDM) to match the normalization on large scales as determined 
from the COBE measurements of the Cosmic Microwave Background (Bunn & White 1997), and a 
flat critical density mixed dark matter model with a hot component added to the CDM (CHDM). 
There are several important and well-established astrophysical measurements which constrain the 
various combinations of cosmological parameters. The parameters for each model, which we list in 
Table ||, have been determined to provide good consistency with these observations. For example, 
the combination Q.sh'^ is restricted by Big Bang nucleosynthesis constraints and the measured 
abundance of primordial deuterium to lie in the range 0.015-0.025 (Copi, Schramm, & Turner 
1995; Buries &: Tytler 1998). In addition, because the H i column density scales approximately as 
(Oft/i^)^ for a fixed UV radiation intensity, we choose O^ and h so that Q,i,h? is the same for three 
of the models sCDM, ACDM and OCDM. The fluctuation normalization in a sphere of 8/i^^ Mpc 
is defined to match observations of the number density of galaxy clusters (White, Efstathiou, & 
Frenk 1993; Bond & Myers 1996) in all the models. In addition, a tilt has been applied to the 
CDM power spectrum in the tCDM model in order to approximately match the amplitude of the 
CMB quadrupole as measured by COBE (Bunn &: White 1997). The cosmological constant in 
the ACDM case is consistent with the upper limit (Oa < 0.7) of Maoz and Rix (1993) and the 
best fit parameters of Ostriker and Steinhardt (1995). One of the major problems with the sCDM 
model is its difficulty in matching observations of the large scale structures in the universe. Since 
the standard CDM model is historically one of the most studied models, however, we use it as 
our canonical model to which the perhaps more viable additional models considered here may be 
compared and through which we investigate the dependence of the Lya statistics on the fiuctuation 
power spectrum. We refer the reader to Zhang et al. (1995, 1997, 1998) for further details and 
results from our previous sCDM simulations. 

The initial data were generated using COSMICS (Bertschinger 1995) with the BBKS transfer 
function (Bardeen et al. 1986) to compute the starting redshifts and the initial particle positions 
and velocity perturbations appropriate for all models except CHDM. We used CMBFAST (Seljak 
&: Zaldarriaga 1996) to solve the linearized Boltzmann equations to set the initial conditions for 
CHDM. For the comoving box size adopted (9.6 Mpc) and the corresponding comoving grid cell size 
(37.5 kpc in our high resolution runs), the relevant wavenumber domain of the simulations at z = 



is 168 > k > 0.65 Mpc^^, where k = 2Tr/i and i is the length scale. Over this domain, the sCDM, 
ACDM and OCDM models all have a similar power distribution. The tCDM and CHDM models, 
on the other hand, have an overall lower normalization (see Table |^) in addition to a steeper slope 
that drops slightly more sharply than the other models over the smaller scales. In Figure ^ we 
show the linear power spectra for these models evolved to z = 3, the redshift at which we present 
many of our results. Since previous work (Zhang et al. 1998) indicates that the sizes of the low 
column density absorbers at z ~ 3 are ~ 100 kpc, it is useful to characterize the models in terms of 
their power at these small fluctuation scales. A useful measure of this power, introduced by Gnedin 
(1998), is 

al, = fpik)e-''''/'l^^ (1) 

(where k = 27r/£, P{k) is the linear power spectrum at z = 3 and ks4 = 34J7g h Mpc~^). This is 
also listed for each model in Table |l|. 

In addition to specifying a cosmological model, it is also necessary to include a background 
UV radiation field to ionize the IGM. We have implemented the spectrum computed by Haardt & 
Madau (1996) for a flat universe on the basis of radiation transfer in a clumpy universe and the 
measured luminosity function of QSOs, accounting for QSO sources, absorption by the Lya forest 
and Lyman limit systems, and re-emission of the recombination radiation from the absorbing clouds. 
This spectrum reionizes the universe between redshifts 7 and 6 and peaks at about z = 2. Because 
the OCDM model corresponds to a cosmology not considered by Haardt &: Madau (1996), we use 
the field for a flat universe for this model also, noting that the neutral fractions are in any case 
rescaled, as described below, to match observations. We also note that only clouds which are 
optically thin at the Lyman edge are considered in this paper. Hence the optically thin limit is a 
good approximation and it is not necessary to account for self-shielding and radiative transfer of 
the external ionizing radiation field. 

The numerical computations were performed using two different numerical codes, Kronos and 
Hercules, each in a simulation box of length 9.6 Mpc comoving with the universal expansion. Kronos 
(Bryan et al. 1995) is a single grid Eulerian code that uses a particle-mesh algorithm to follow the 
dark matter and the piecewise parabolic method (PPM) to simulate the gas dynamics. Since non- 
equilibrium chemistry and cooling processes can be important, six particle species (H i , H ii , He i , 
He II , He iii and the electron density) are followed with a sub-stepped backward finite-difference 
technique (Abel et al. 1997; Anninos et al. 1997). This is the same non-equilibrium chemistry 
and cooling model used in our previous studies of the Lya Forest (Zhang et al. 1995, 1997, 1998; 
Charlton et al. 1997; Bryan et al. 1998). For sCDM, ACDM , OCDM, and tCDM we use 256^ grid 
cells in the simulation box to follow the evolution of 128^ dark matter particles. This results in a 
spatial resolution of Ax = 37.5 kpc. For CHDM 128'^ grid cells are used with 64^(128^) cold (hot) 
dark matter particles, respectively, resulting in a lower spatial resolution. Ax = 75 kpc, for this 
model. For the sCDM simulations with varying (Tg/j-i we present simulations with both 128'^ and 
256^ grid cells resulting in both low and high spatial resolutions of 75 kpc and 37.5 kpc, respectively. 



We have also simulated three of our models (sCDM, ACDM , OCDM) with a different numerical 
code, Hercules (Anninos, Norman & Clarke 1994; Anninos et al. 1997). Hercules is a nested grid 
code that utilizes a multiscale PM method for the dark matter, artificial viscosity methods for 
the baryonic fluid, and the same non-equilibrium chemistry and cooling model as above. The 
simulations produced from this code use 128^ particles and 128'^ cells for both the nested and 
parent grids. However, in order to derive a more representative sample for statistics, the results 
discussed in this paper are extracted from the parent grid only. Thus these simulations are of lower 
spatial resolution than most of the Kronos simulations, although the dark matter mass resolution 
is the same. For statistics that are insensitive to spatial resolution, a comparison of the results of 
the two codes is useful to insure that simulation results are robust against changes in numerical 
technique. 

Synthetic spectra are generated along 300(900) random lines of sight through the Kronos 
(Hercules) simulated volume using the method of Zhang et al. (1997) including the effects of 
peculiar velocity and thermal broadening of the gas. (We have verified that decreasing the sample 
size from 900 to 300 for the Kronos data does not affect the results except for a slight increase in the 
scatter of the line properties for the highest, optically thick column density systems, a regime where 
our results become unreliable anyway because of the absence of radiative transfer in the codes.) 
Since we are primarily concerned in this paper with a comparison of model predictions, we have 
not included noise or continuum fitting in the analysis. Furthermore the resolution of the spectra, 
1.2km/s, is the same for all the simulations, a value that is smaller than current observations. 
However, we have shown elsewhere (Bryan et al. 1998) that, as long as we restrict ourselves to 
high quality observational data, the impact of not including these observational difficulties is small. 
In addition to analyzing the raw optical depth and flux distributions, line lists are extracted from 
the data using a Voigt profile fitting procedure. This is described in more detail elsewhere (Zhang 
et al. 1997), but we outline it briefly here. First, maxima in the optical depth distribution are 
identified as line centers. Then Voigt profiles are fit, using a non-linear minimization, to the part 
of the spectrum which is above thi = 0.05 and between neighbouring minima. This results in the 
same spectral threshold Ft = e""^* = 0.95 as the high resolution Keck HIRES spectrometer. Each 
line of sight chosen produces a sample spectrum with on the order of 10-100 lines per redshift 
interval (5z = 0.1 depending on the redshift and cosmological model. The statistics of these linelists 
are discussed in §^ and §|5|. 

The amplitudes of the distributions found in the models cannot be used as a basis for comparing 
the models since they may be arbitrarily re-scaled for any individual model using the ionization bias 
factor b\oa = ^^/F, where Q-b is the fraction of the critical density carried in baryons and F is the 
Haardt-Madau (1996) parameterization of the metagalactic UV ionizing background extracted from 
the observed distribution of quasars. It is important to normalize all the models consistently before 
comparing the shapes of any of the distributions. This may be done in a variety of ways. We do so 
by matching the mean H I opacity in each simulation to the measured intergalactic H I opacity at 
2 = 3. In Zhang et al. (1997), we found the opacity measurements of Steidel & Sargent (1987) and 



Zuo & Lin (1993) gave a mean H i opacity at z = 3 of Tq, = 0.27 — 0.35, although values as much as 
30-60% larger have been claimed (Press et al. 1993; Ranch et al. 1997). Because of the uncertainty 
in this measurement, we also require consistency with the number density of lines observed above 
a threshold of log A'hi = 13.5. Using the three quasars in Hu et al. (1995) for which lines in the 
full redshift range 3 < z < 3.1 are listed, we find a total of 61 lines for the three lines-of-sight in 
this redshift interval with log A^hi > 13.5, for which the line lists should be complete. (An estimate 
based on using the available lines for all four QSO line lists in Hu et al. in the redshift interval 
2.9 < z < 3.1 gave essentially the same line density.) Normalized to fa = 0.30, the CHDM, sCDM, 
ACDM , OCDM, and tCDM models predict, respectively, 60.8, 62.1, 62.7, 63.7, and 59.5 lines, in 
close agreement with the observed number. Normalizing to fa = 0.35, the respective numbers of 
predicted lines are 73.7, 74.3, 73.9, 75.2, and 72.8. While these are not badly inconsistent with 
the observed number, they are all fairly high. We normalize the spectra according to Tq, = 0.30 
throughout this paper, noting that this value is still not well agreed upon. In Figure ^ we plot a 
related statistic, Te// (Zhang et al. 1997) for the normalized spectra of our models and compare 
to recent data by Kirkman & Tytler (1997). After normalization all of our models are consistent 
with the data over the redshift range 2 < z < 4 considered by this paper. 



3. Direct Optical Depth and Flux Measurements 

Historically Lya absorption spectra have been analyzed in terms of the statistics of spectral 
line features and as such have been plagued with difficulties of the line fitting procedure such as line 
identification and blending. Many of these difficulties become increasingly severe at higher redshifts 
making the results of the analysis uncertain. It is thus natural to ask whether statistics dependent 
directly on the observed flux and optical depth without recourse to line fitting might be of use 
in describing the forest and discriminating among competing models. Statistics of this kind have 
recently been proposed by several authors (Miralda-Escude et al. 1997; Ranch et al. 1997; Cen 
1997). Since these nonparametric measures are also easier to relate theoretically to the physical 
state of the absorbing gas, we begin our discussion with them. 



3.1. Optical Depth Probability Distribution Function 

The optical depth r is related to the transmitted flux F hy F = exp{—T). We define the optical 
depth probability distribution dP/dr as the probability that a pixel will have optical depth between 
T and r + dr. In Figure ^ we use spectra generated from the sCDM high resolution simulation to 
show rdP/dr versus r for redshifts 2 = 2,3 and 4 (top panel). Although the peak of the distribution 
decreases and the distribution broadens slightly with decreasing redshift, the principal contributor 
to the redshift evolution seen in Figure |3| is the evolution of the optical depth r. Hui, Gnedin & 
Zhang (1997) discuss in detail the dependence of r on the distribution and properties of neutral 
hydrogen along the line of sight in an expanding universe. Since we would like to understand the 



redshift evolution of the optical depth in terms of simple scaling laws, we repeat some of their 
discussion here in order to the isolate the key factors controlling this redshift evolution and clarify 
the scaling law assumptions. The optical depth is defined as 

. . f^b dx . . 

t{uo) = / nniCTa r—— (2) 

Jxa i + 2 

where Vq is the observed frequency, nni is the number density of neutral hydrogen, z is the redshift 
of the absorbing gas, cJq, is the absorption cross section for Lya , and the integral is over the line of 
sight between the quasar (xa) and the observer (xb) in comoving coordinates. In practice the form 
of the Lya absorption cross section limits the integration range per absorber to a small portion of 
the line of sight. It is thus useful to make a change of variable to velocity coordinates u about some 
characteristic average redshift z in the problem. For example, for simulated data the redshift z 
might be a given output redshift for the simulation. The observed frequency Uq and the frequency 
u of the radiation in the absorber rest frame are then related by 

u = Uoil + z){l + u/c) (3) 

where 

H(z)(x - x) . . ,^, 

U = "■ '"■ + Vpec{x). (4) 

1 I z 

X is the comoving position along the line of sight whose redshift is exactly z, Vpec is the physical 
velocity of the gas, and H{z) is the Hubble parameter defined by 



H{z) = Ho^jQmil + zf + {l-nm- Oa)(1 + zy + Ua. (5) 

The first term in Equation |^ represents the contribution of the residual Hubble flow about the 
mean while the second term is due to the physical bulk flow of the gas. We assume u/c « 1 and 
neglect contributions from turbulent flows since they would be unlikely in the low column density 
regions we are considering. Under this change of variable the Lya cross section becomes 

where aao = 4.5 x 10^^*^ cm^ sets the scale of the absorption cross section in terms of fundamental 
constants, uq is the velocity u for which the frequency i^ in the rest frame of the absorbing gas is 
equal to the Lya frequency v^ and b = J2kBT/mp is the thermal width. For absorption lines of 
neutral hydrogen with column densities A'hi < 10^^ cm~^ the thermal profile dominates the cross 
section so we neglect the contribution of the natural line width to a a ■ The optical depth r can now 
be written as 

^-{u-uoY/h' dn (7) 



streams 



si.rfam.s V ^^ / 



du 



dx 



The sum over streams represents the possibility that a given velocity u corresponds to more than 
one position x. Although the integration formally runs over the full line of sight from quasar to 



observer, the Gaussian form for the cross section effectively hmits the u integration to a narrow 
range around uq (thus justifying our replacement of z everywhere by z.) To simplify notation we 
drop the bar letting z represent z in what follows. 

We assume that the number density of hydrogen traces the baryon gas density well. (There 
has been little metal production in these low density regions and there is no interaction that would 
cause the helium and hydrogen to separate.) Thus the number density of neutral hydrogen is 
riHi = PhXni where Xhi is the neutral fraction and p^ the gas density. In ionization equilibrium 
(which is well satisfied except for the period of initial reionization) the neutral fraction of hydrogen 
is Xhi oc p^T^^ such that the number density of neutral hydrogen (relevant to Lya absorption) 
scales as 

nni oc {nyfT-Hz){l + zf{l + 6,fT-^-' (8) 

where (5f, is the baryon overdensity. Studies (Hui & Gnedin, 1996; Weinberg et al. , 1996) of the 
equation of state for the gas find that for unshocked gas at low to moderate baryon overdensities 
{^b ^ 5) the equation of state is well fit by a power law: 

T^{l + zf-\l + 6,y-\ (9) 

Thus 

UHi oc {n,h^fT-\z){l + zf-\l + (^,)2-7-0-7^. (10) 

For a uniform radiation field and reionization that occurs before z = 5, as is the case in our 
simulations, 7 ~ 1.4. This is in agreement with the value 7 ~ 1.5 found by Zhang et al. (1998) 
for clouds with column densities in the range 12.5 < log A^hi < 14.5. Furthermore the assumption 
that most of the optical depth arises from low column density absorbers, large structures whose 
overdensities and peculiar velocities are slowly varying compared to the thermal profiles, means that 
multiple streaming is rare, the sum over streams in Equation ^ can be dropped, and 
We then integrate over the thermal profile to obtain (Croft et al. 1997) 



H 

1 + 2- 



ro.'^^^^§^{l + zr{l + 5,r. (11) 

Note that in this limit r need no longer have a thermal profile about its maximum (Hui, Gnedin 
& Zhang, 1997). If (5^ is evolving slowly over this redshift range, r should scale as 

In the middle panel of Figure |^ we use this simple scaling law to rescale the z = 2 and z = A sCDM 
distributions from the top panel to z = 3, the redshift at which all the models are normalized. We 
do this in order to test how well the simulations obey this simple scaling relation: if they followed 
it exactly then all three curves would overlap. Most, but not all, of the redshift evolution of this 



distribution is accounted for by the scaling of r given in Equation 12. Since the evolution of the 



metagalactic UV radiation field F is relatively slight over this redshift range, we are left with the 



remarkable conclusion that most of the evolution of the Lya forest is a direct consequence of the 
universal expansion. The direct numerical results of Zhang et al. (1998) support this conclusion. If 
we include the evolution of the baryon overdensity, as shown in Figure ^ for the sCDM simulation, 
and shift the overdensity distribution until the peaks overlap, the (l+^b)^''^ dependence in Equation 
pi] for the optical depth distributions predicts additional scaling factors of ~ 1.64 (0.77) at z = 2 
{z = 4), respectively, for r that bring the distributions (shown in the bottom panel of Figure |3|) 
into close agreement. The remaining small differences, the slight broadening of the distribution 
and a reduction in its peak amplitude with decreasing redshift, most probably reflect the fact that 
the shape of the baryon overdensity distribution is also evolving slowly with z. 

The top panel of Figure |^, shows rdP/dr for the five sCDM simulations with varying power 
and spatial resolution. From this we can see that the optical depth PDF at a given redshift is 
insensitive to the spatial resolution of the simulation. The shape of the distribution, however, 
is strongly dependent on the amount of small scale power present. Models with less power at 
these scales produce narrow, sharply peaked distributions. As the power increases, the distribution 
flattens and broadens. In the lower panel of Figure || we show rdP/dr versus r for the simulations 
in the model comparison study. These distributions again display a clear dependence on the power 
spectrum of the model with OCDM (our model with the most small-scale power) producing the 
broadest distribution, and CHDM and tCDM (models with the least small-scale power) producing 
the most sharply peaked distributions. Thus this statistic is particularly promising as a model 
discriminator in that these differences between models are significant in the range 0.02 < r < 4 
that should be accessible to observers. 

We quantify this relation between the shape of the r distribution and the amplitude of the 
power spectrum by fitting a log-normal to the curves: 

dr 

Although this does not fit the profiles in Figure ^ in detail, it does provide an adequate description 
as long as we restrict the range of optical depths used in the fitting. Here we adopt 0.02 < r < 4, 
corresponding roughly to the observable range. A different range or a different fitting function 
changes the details, but not the nature of our result. In Figure ^, we show the correlation between 
(Tt, a measure of the width of the distribution, and (T34, the amplitude of the linear power spectrum 
on small scales as defined in Equation |l|. The strength of the correlation is striking. The low 
scatter around the power law relation shown in this figure bolsters our claim that the shape of 
the r distribution function is insensitive to other cosmological parameters. To give an idea of 
the uncertainty in each point, we fit both the high and low resolution simulations for the sCDM 
o"8 = 0.3 and as = 0.7 models. In both of these cases o",- differs by less than 10 %. 
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3.2. Flux Probability Distribution 

Although the optical depth PDF is easier to model theoretically, the flux PDF (where dP/dF 
is the probability that a pixel will have transmitted flux between F and F + dF) is closer to actual 
observation. The top panel of Figure shows the flux probability distribution functions for the 
high spatial resolution sCDM model with erg = 0.7 at z = 2, 3 and 4. The bottom panel of Figure 
shows the prediction of the simple r scaling given in Equation 12 applied to the flux and these 



same flux probability distributions. Again we attempt to rescale the z = A and z = 2 distributions 
to z = 3 in order to test the scaling. This results in a highly nonlinear mapping of the flux and the 
flux PDF from z = j to z = ^ given by dPj/dF -^ TjF'^-'^/'^dPj/dF and Fj -^ FJ^'' for r,- -^ tj/tj, 
where t] = 0.356(6.511) for j = 2(4), respectively. While the shapes of the distributions in the 
top panel appear quite different, much of the z evolution of the flux probability distributions is 
explained by this simple scaling, the remainder representing mostly the effect of the evolution of 
the baryon density in the cosmological model. We do not plot the scaled distribution for z = 4 
below the scaled flux of 0.5 because this already corresponds to an unsealed flux of 0.015, close to 
saturation and most likely noise dominated in the observations. 

The flux PDF depends only weakly on simulation grid resolution (Bryan et al. 1998). Its shape 
is strongly dependent on the power spectrum of the underlying cosmology. In the top panel of Figure 
^ we show the flux probability distributions in the sCDM model (spatial resolution Ax = 37.5 kpc) 
for cluster scale normalizations cxg = 0.3 and 0.7. The dependence on the normalization of the power 
spectrum is clear. The number of pixels found with flux in the central flux range 0.3 < F < 0.9 
is greater for models with less power (ag = 0.3); while the number of pixels with flux in the low 
{F < 0.3) and high (F > 0.9) flux ends of the distribution are less than for models with greater 
power (cTg = 0.7). This is in qualitative agreement with Croft et al. (1997a). We note, however, 
that our result (using the Kronos code) does not require any smoothing of the simulations as was 
the case for their TreeSPH simulations. In the lower panel of Figure |8| we present the z = 3 flux 
PDFs for the five models of the model comparison study. Models with lower power at small scales 
(tCDM, CHDM) have a larger flux PDF for 0.3 < F < 0.9 than sCDM and ACDM , while the low 
density model (OCDM) with the highest spectral power at these scales has the smallest flux PDF 
in this range, as expected. Furthermore, the differences between models can be substantial. For 
example, at F = 0.6 the OCDM results he 10 % below the ACDM model result while the CHDM 
result lies above the ACDM result by about a factor of 1.4. We remind the reader that the mean 
of the distribution has been fixed to match observations. Thus this statistic should be useful to 
constrain competing models. 



3.3. Fraction of high Lya opacity 

Another possibly useful statistic for discriminating models is the fraction of a quasar spectrum 
with Lya optical depth greater than a specified value tq, i.e. the cumulative distribution in optical 
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depth. Small differences in the amplitude normalization of the primordial power spectrum may be 
enhanced in the cumulative opacity data (Cen 1997). 

Figure ^ shows the linear correlation between the opacity at line center and the column density 
of absorption features in the sCDM model, ranging from the optically thin to thick at the Lyman 
edge. The nearly unbroken relation Tc oc Nhi, which exists down to the incompleteness density 
of ~ 10^^ cm~^ is attributed to the weak correlation between the Doppler parameter and column 
density since, in general, Nhi oc bTc- The lower bound on the opacity {tc > 0.05) is set by the 
transmission or spectral threshold Ft = e""^* = 0.95 used in the line identification procedure. Using 
Figure ^ as a guide, we investigate the cumulative opacity distribution with the following minimum 
opacity thresholds: r > 0.1, 1, and 7 which, if associated with the line centers, would correspond 
roughly to column densities of log Nhi = 12.5, 13.5 and 14.5 respectively. The distributions P{t > 
tq) for the above minimum opacity thresholds are plotted in Figure |l^ at redshifts z = 2, 3, and 4 
for the models in the model comparison study. In comparing groups with the same tq, the smaller 
threshold curves are more highly clustered and less sensitive to the background cosmological model 
parameters. This is especially evident in Figure ^ where we show the cumulative distributions of 
the optical depth at redshift z = 3 for these models. 



4. Line Parameter Statistics 

In this section we present a line analysis of the spectra generated from the various model 
simulations. We compare and contrast the cosmological models based on the column density 
distribution and the evolution of line number. 



4.1. Hi Column Density Distribution 

One of the most robust line statistics used in the analysis of the Lya forest is the H i column 
density distribution, which is well converged by simulation box sizes of 9.6 Mpc and is insensitive 
to changes in the simulation grid resolution or treatment of gas hydrodynamics (Bryan et al. 1998, 
Zhang et al. 1997). The H i column density distribution, defined to be Nhi = Jx^ j?^dx, is closely 
related to the optical depth r through the dependence of each on the number density of neutral 



hydrogen. Thus using Equation 10 and the approximations that led to Equation |ll| we expect the 
H I column density to scale as 



Nhi oc 



T{z)H 



{l + zt^ j{l + 5bf-Uu (14) 



In the top panel of Figure ^ we show the raw (uncut) H i column density distribution for the high 
resolution sCDM model for redshifts z = 2,2> and 4. In the bottom panel of Figure ^ we see that 
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the column density scaling 

iV/./ocip^(l + .r, (15) 



the same relation as the naive scaling relation given in Equation 12 for the optical depth r, accounts 
for the redshift evolution of the column density distribution amazingly well. This demonstrates 
that the column density, an integrated quantity, is much less sensitive than the optical depth 
distributions to the redshift evolution of the gas overdensity within an absorbing structure. The 
differences seen in the low column density end of the distributions, particularly for z = 4, may be a 
result of the simulation spatial resolution (Bryan et al. 1998), while the differences observed in the 
high column density end may partially be due to shot noise in the high z data. In Figure O (top) 
we explore the dependence of this distribution at a given redshift (z = 3) on the power spectrum 
of the underlying cosmology. We find qualitative agreement with semi-analytic arguments (Gnedin 
1998; Hui, Gnedin, & Zhang 1997) in that models with less power on small scales (such as sCDM 
with a% = 0.3 and C734 = 0.812) have H i column distributions with significantly steeper slopes than 
models (such as sCDM with as = 0.7 and (T34 = 1.89) with more power at these scales. However, 
as we discuss in more detail below and in Table ^, quantitative agreement between the simulations 
and the predictions of these semi-analytic arguments seems more difficult to achieve. 

In Figure ^ (bottom), we show the H i column density distribution at redshift z = 3 for 
(Kronos) simulated spectra in the model comparison study and compare the simulated data with 
data from Kirkman & Tytler (1997) and the fits provided by Kim et al. (1997). The distributions 
are conventionally quantified by fitting them to power laws, dN/dNm oc N^-[ . We use the same 
sets of column density cuts on the simulated data in the model comparison study as Kim et al. in 
order to expedite comparison with the data and use a direct unweighted least squares fit (all quoted 
errors are 2a) to extract the slope /3 from the simulated data. Our results are summarized in Table 
y. We find again the expected dependence on the fluctuation power spectrum. For the column 
density range 13.7 < log A^hi < 14.3 (given by the column labeled /3h) the most shallow slope is for 
OCDM, the low matter density model with 0-34 = 2.50 while CHDM and tCDM with 0-34 = 1.14 
and 1.09, respectively, give the steepest distributions (see Table |l|). The predicted column density 
distributions generally also steepen with time (decreasing redshift). Kim et al. find /3 = 1.46±0.07 
(2a") for this column density range at {z) = 2.85. This is formally inconsistent with all the models 
at the 3a level except for OCDM, although it is marginally consistent with ACDM . 

Results for the column density range 12.8 < log A'^hi < 14.3 in Table are shown in the 
column labeled Pi. The average distributions are generally shallower when extended to lower 
column densities, showing that the distributions are curved. Kim et al. similarly find a shallower 
distribution over this column density range with their results for lines at {z) = 2.31, {z) = 2.85, 
and (z) = 3.35 shown in the last row of Table |2[ The result at (z) = 2.31 is inconsistent with 
all of the simulation results at z = 2, but note that the quoted uncertainty in the observation is 
eight times smaller than at (z) = 2.85, despite comparable numbers of absorbers. By contrast, 
the result at (z) = 2.85 is formally consistent at the 3<t level with all the models. At (z) = 3.35, 
Kim et al. find /? = 1.59 ib 0.13 {2a). This result is consistent with the simulation results for 
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the OCDM model and marginally consistent (at the 3cr level) for the sCDM and ACDM models. 
The observational data, however, also suggest a weak steepening of the distribution with increasing 
redshift, contrary to our findings. These discrepancies might indicate that the redshift evolution of 
the ionizing radiation field may be somewhat different from that of the Haardt & Madau spectrum 
assumed in the simulations. 

We may compare the simulation results with the semi-analytic predictions of Hui, Gnedin, &: 
Zhang (1997) to understand the trend of changing steepness with power spectrum. We provide the 
predicted values of /3 according to the prescription of Hui et al. in Table ||. We assume T oc p^^, as 
found by Zhang et al. (1998) for this column density range in an sCDM simulation. The uncertainty 
in the T — pB relation introduces only a 10% uncertainty in the prediction for /3 — 1, so it seems 
reasonable to retain it for the other models as well for this purpose. The predicted values of (3 
for the sCDM, tCDM, and CHDM models at z = 3 match the simulation values to within la, in 
agreement with the comparison in Hui et al. with one of our earlier sCDM models. However, the 
predictions for OCDM and ACDM at z = 3 are in disagreement with the semi-analytic arguments 
giving too steep a slope. The predictions do less well for all models at z = 2. In particular, the 
simulation results show a steepening of the column density distribution toward decreasing redshifts 
for all the models, opposite to the predicted trend. 

Over the wider column density range 10^^'^ < A'hi < 10^^ (summarized in the column labeled 
/5f in Table |2|), we see that the average distributions continue to steepen toward higher column 
densities. Kim et al. obtain (3 = 1.46 for this column density range at {z) = 2.85, with a steepening 
to /? = 1.55 at {z) = 3.7. The results for the tCDM and CHDM models (1.95 ±0.06 and 1.92 ±0.06 
at z = 3, respectively) are substantially steeper than these values. Because the distribution deviates 
from a pure power law at the low column density end, it is useful to split the simulation samples 
into two halves, fitting each to a power law. These results are given in the last two columns of 
Table ^ where /3i is the slope of the column density distribution for lines at the low column density 
end (lO^^*^ < A^Hi < 10^^) and /32 is the slope of the column density distribution for lines at the 
high end (10^^ < iVni < 10^^). Giallongo et al. (1996) obtain /3 = 1.8 for systems with iVm > 10^^ 
and 2.8 < z < 4.1. 

Finally we note that the analogous column density distributions derived from the lower reso- 
lution Hercules runs give similar results and slopes as the Kronos data. For example over the full 
column density range at z = 3, Hercules data give slopes for the column density distribution of 
1.71, 1.66 and 1.62 for the sCDM, ACDM and OCDM models, respectively, consistent within errors 
with the Kronos results. This suggests that the distribution function is a robust diagnostic, being 
relatively insensitive to grid resolution and numerical method. A preliminary comparison with 
the data favors models with more power at these scales than in our CHDM or tCDM cosmologies. 
However, there appears to be some discordance in the observations, so a more definitive comparison 
will require more work. 
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4.2. Line Number Evolution 

The number of Lya lines at a particular redshift reveals how many intergalactic absorbers exist 
at that time between the quasar and observer and, given certain assumptions on their geometry, 
the size and volume filling factor of the absorbers can also be deduced. Since the column density 
of Lya lines corresponds to the mean overdensity and size of the clouds fairly well (Charlton 
et al. 1997; Zhang et al. 1997), it is useful to see how the number of lines evolves with different 
column density cutoffs, as this will track the evolution of morphologically distinct small scale 
structures in the universe. 



Figure 14 shows the evolution of the number of lines with H i column densities greater than 10^^ 
cm~^ , 10^^'^ cm~^ , and lO^'* cm~^ , respectively, comparing results for the models in the model 
comparison study with the observed data from Kulkarni et al. (1996) at z ~ 2, Hu et al. (1995) 
at z ~ 3 and Lu et al. (1997) at z ~ 4. For a fixed transmission cutoff (here Ft = 0.95) and 
column density threshold, the total number of lines per unit redshift decreases with time because 
the opacity of the universe decreases from both the increasing flux of radiation and the expansion of 
the universe. With the exceptions of the tCDM and CHDM models for the lowest column density 
threshold A^hi > 10^^ where incompleteness due to line blending becomes significant at higher z, 
the deviation from a fixed power law behavior tracks predominantly the behavior of the radiation 
flux. Fitting the evolutions to the form dN/dz oc (1 + z)'^ over the range 2 < z < 4, we find the 
exponents are fairly similar in the different models. We summarize these results in Table |3| (all 
errors are 2a). To compare these simulated results with observational data we fit the combined 
line lists from Kulkarni et al. , Hu et al. , Kirkman &: Tytler, and Lu et al. to the same power 
law behavior and display those results in the row labeled "combined" in Table ^ again with 2a 
errors. (Lines near the QSO emission redshift were avoided because of the proximity effect, as were 
lines associated with metal systems.) Kim et al. obtain 7 = 2.78 it 1.42 {2a), fit over 2 < z < 3.5 
for systems with column densities lO^'^'^^ < A^hi < 10^^ cm~^ . Using the same column density 
cuts and simulation data for z = 2 and 3 only, we find power law exponents (labeled 7h in Table 
I) for the sCDM, ACDM , OCDM, tCDM and CHDM models of the comparison study in good 
agreement with the observational result. For the lower column density range lO^'^'^ < A^hi < 10^^ 
cm~^ , the power law exponents are labeled 7^ in Table y. Kim et al. (1997) obtain 7 = 1.29 it 0.90 
{2a) fit over 2 < 2; < 4. Our sCDM, ACDM and OCDM model predictions are in good agreement 
with this observational result, although tCDM and CHDM show a somewhat stronger evolution. 

All the models yield comparable levels of evolution at each of the column density cutoffs and, in 
fact, the evolution slopes in all the models agree for the most part within errors with the observed 
values. Two trends in the model predictions are apparent. The first is that at a given column 
density threshold the slope 7 of the line number evolution is correlated to the slope /5 of the column 
density distribution with 7 increasing for models with larger (5 (i.e. models with less power in the 
fiuctuation spectrum at small scales). The second trend is that stronger column density lines for a 
given model exhibit a greater rate of redshift evolution. This is not unexpected since the evolution 
for a fixed transmission threshold is essentially determined by the radiation field (which is the same 
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in all the models) and, to a lesser degree, the expansion of the universe. In our previous studies 
(Zhang et al. 1997) we have found little intrinsic cloud evolution over these redshift intervals. If we 
assume that the evolution of the gas overdensity does not contribute significantly to the evolution 
of line number or column density at these redshifts we can use Equation |l^ and the power law 
dependence of the column density distribution dN / [dN-mdz) oc -/Vjjj to predict the number of lines 
above a fixed column density threshold as a function of redshift z. We find 

dn ( (1 + zY-^V^'^'^ 



In Figure |15| we compare the scaling predictions of Equation |16| to the simulation results using the 
high resolution sCDM erg = 0.7 model for column density thresholds A'^hi > 10^^, 10^^'^, and 10^^ 
cm~^ , respectively. Since /3 also evolves weakly with z, we use /? at z = 3 as representative of the 
average for 2 < z < 4 in Equation |l6[ For the lowest two column density thresholds we use the single 
power law fit to the column density distribution over the range 10^^'* < A^hi < 10^^ cm~^ , while 
we use (32 from the two power law fit for the high (A^hi > 10^^ cm~^ ) column density threshold. We 
normalize the scaling predictions to the simulated number of lines at z = 3 because that is where 
all the models in our study were normalized to the observational data. For the lower two column 
density thresholds, the scaling prediction tends to overestimate the number of lines at z = 4. For 
the lowest column density threshold (A^hi > 10^^) this again is partly due to incompleteness in the 
simulation line lists caused by line blending, an effect that becomes more severe for low column 
densities at high z. Furthermore the slope of the low column density end softens for A'^hi < 10^^ 
with the break at A'hi ~ 10^^ probably reflecting a change in the absorbers from low density 
structures evolving primarily with the universal expansion to structures undergoing gravitational 
collapse (Bryan et al. 1998). This deviation of the column density distribution from the pure 
power law assumed in the scaling relation would also cause the scaling law to overproduce the low 
column density lines. For the column density threshold A'hi > 10^^'^ the discrepancy between the 
scaling prediction and the simulation results at z = 4 is reduced. This is to be expected since the 
high spatial resolution sCDM linelists should be complete at this column density threshold. For 
lines with A^hi > 10^^, where the absorbers share a common morphological type and the column 
density distribution is well fit by a single power law, agreement between the scaling prediction and 
the simulations is good. In the lower panel of Figure |l^ we compare the scaling predictions with 
the simulation results for the models in the model comparison study for the high column density 
threshold case. The scaling predictions for all models agree reasonably well with the simulations. 



Encouraged by these results we solve Equation 16 for the shape of the UV ionizing background 
r(z) in terms of the Hubble parameter H{z) modeling the universal expansion and the (in principle) 
measurable quantities dn/dz{Niii > W^'^) and P2, the slope of the column density distribution over 
this column density range. 



We use Equation 17 to compute r(2;) with simulation data from the model comparison study and 
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compare these predictions to the Haardt-Madau spectrum actually used in Figure |T6|. Although 
the prediction is highly sensitive to the slope of the column density distribution used (whose errors 
are still quite large), it is gratifying, given the simplicity of the scaling relations, that all of the 
models reproduce the assumed Haardt-Madau evolutionary trend. 



5. Doppler b Parameter 

Recent papers (Bryan et al. 1998, Theuns et al. 1998) have shown that both the Doppler b 
parameter and a related nonparametric statistic, the mean flux difference as a function of velocity, 
require high simulation spatial resolution to model properly. In this section we investigate the 
dependencies of these statistics on the properties of the cosmological model. Because these statistics 
are highly sensitive to the spatial resolution of the simulation, we present results only for those 
models (sCDM, ACDM , OCDM, and tCDM) simulated with our highest spatial resolution Ax = 
37.5 kpc. 

The Doppler b parameter measures the amount of line broadening due to thermal broaden- 
ing, physical velocities, Hubble expansion broadening and the shape of the absorber density profile 
(Bryan et al. 1998). Both Hubble and thermal broadening are significant for the lower column 
density lines that arise from structures found in voids that are still expanding in absolute coordi- 
nates. The thermal contribution only becomes dominant for the higher column density lines that 
have turned around and are gravitationally collapsing. Furthermore, the b parameter is highly 
sensitive to the simulation spatial resolution. Lower resolution simulations numerically thicken the 
lines causing the width of the lines, and thus b to be overestimated (Bryan et al. 1998; Theuns 
et al. 1998). In our previous work (Bryan et al. 1998) we argued that the shape of the 6-distribution 
was in rough agreement with observation and particularly that the high b power law tail of the dis- 
tributions arises naturally in hierarchical models when quasar lines of sight pass obliquely through 
the filamentary absorbing structures (Rutledge 1998). However the median of the simulated b pa- 
rameter distribution for the sCDM model, calculated from simulations with high spatial resolution, 
was now substantially below the ~ 30 km/s median seen in the observations. Thus the sCDM 
model, which previously had appeared to be in agreement with observations, is now discrepant. In 
Figure |17| we show the Doppler b parameter distributions extracted from the high grid resolution 
(Ax = 37.5 kpc) Kronos simulations at redshift z = 3 for the sCDM, ACDM , tCDM, and OCDM 
models for lines with column densities between 10^^'^ cm~^ < A'^hi < 10^^ cm~^ . We present for 
comparison data from Kim et al. (1997) for z = 3.35. The ACDM and OCDM models, like sCDM, 
have their b distributions shifted too much to the left (to low b values) to agree with observation. 
Only tCDM, the model with the least fluctuation power at small scales and thus broader density 
structures at this redshift, has a median b approaching the observational values. We explore this 
dependence on the fluctuation power spectrum with the two highest resolution sCDM models (with 



(Tg = 0.7 and 0.3, respectively) in the lower panel of Figure 17 and see that indeed the model with 



lower spectral power produces a b parameter distribution shifted towards higher b (as predicted by 
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Hui & Rutledge 1998). The increase in b for models with less fluctuation power at small scales may 
be partly due to line blending effects at these low column densities. However as shown below, the 
shift to higher b values for these models persists for lines with higher column densities as well where 
line blending should not be as significant and thus can not be explained by line blending alone. 

To facilitate a better comparison of the models with observations we plot the median Doppler 



parameters as a function of redshift in Figure 18 where we have imposed the same column den- 



sity cuts on the lines as those used by Kim et al. (1997). The median b for lines with column 
densities lO^^-*^ cm^^ < A^hi < 10^^ cm^^ and 10^^-^ cm^^ < A^hi < 10^^ cm^^ are shown in 
the top and bottom panels, respectively. While the ACDM , sCDM, and OCDM models predict 
roughly the observed evolutionary trend for both sets of column density cuts, the median b values 
lie systematically more than 6 km/s below the observational data. OCDM, the model with the 
most power at these scales, is the most discrepant. Although tCDM predicts median b parameters 
more consistent with observation, the redshift evolution predicted by this model appears to be in 
disagreement with the data. We can compare these results with other recent data sets. Confin- 
ing to lines with A'^hi > 10^^ cm~^ , we obtain from the published line lists, for 1.9 < z < 2.0 
(1(7 errors), (bmeam ^median) = (32.1 it 2.6,29.7 it 3.3) kms^^ (1(t) (Kulkarni et al. ), 3 < z < 3.1 

(Vean, Vedian) = (38.0± 1.6, 33.6±2.0) kms~l (Hu et al. ), (6„iean,6median) = (27.3± 1.9, 25.9± 1.2) 

kms"^ (Kirkman & Tytler), and at 4 < z < 4.1 (femean, ^median) = (32.6±2.4, 25.9±3.1) kms'^ (Lu 
et al. ), again clearly discrepant with the model predictions. Thus none of the models considered 
here can restore agreement with the observational data. 

We also argued in Bryan et al. (1998) that this discrepancy is not a result of the particular 
choice of line fitting algorithm, but appears for sCDM in the nonparametric moments of the two 
point fiux distribution functions as well. The two-point function P2{Fi,F2, Av) gives the probability 
that two pixels with separation Av will have flux Fi and F2. We plot the normalized moments of 
this function averaged over the flux range Fa to Fi, given by 

Jg dF, J^ dF2P2{F^,F2,Av){F^ - F2) 

j;yFiJ^dF2P2{F,,F2,Av) 

which represents the average flux diffence as a function of velocity for pixels in the range Fa to Fb. 



In Figure 19 we plot the above statistic at 2: = 3 as a function of velocity for several flux ranges for 
the high resolution models (sCDM, ACDM , OCDM, tCDM) of the model comparison study (lower 
panel) and study its dependence on the small scale fluctuation power (top panel) using sCDM with 
as = 0.7 and 0.3, respectively. There is little difference for low velocities, independent of flux level, 
due to the high coherence of the lines. At very large velocity differences there is no coherence and 
the value is just the difference between the mean value of the transmitted flux and the mean flux 
in a given flux interval (Bryan et al. 1998). It is at intermediate velocity separations where the 
statistic is heavily influenced by the structure of the lines. There ACDM , OCDM, and sCDM with 
o"8 = 0.7, whose power spectra are very similar, produce very similar distributions; while the tCDM 
model is quite distinct. We may quantify these model differences by determining at what Av the 
model prediction passes through a given average flux difference. For the flux interval < F < 0.1 
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the simulation predictions pass through the mean flux difference of 0.3 at Av ~ 35 km/s for ACDM 
, OCDM, and sCDM as = 0.7 and at Av w 45 km/s for tCDM. Although observational data is 
limited, these are all lower than the Av ~ 55 km/s from Figure 3 of Miralde-Escude et al. (1997). 

It is important to ask what is needed to restore agreement between the simulations and ob- 
servations. Although we can not completely rule out the possibility that the line fitting algorithm 
contributes to differences in the simulated and observed 6-parameter distributions, we argue that 
its effect should not be significant because the discrepancy is seen at a comparable level (Bryan 
et al. 1998) in the fit-independent two-point distribution of the flux as well. The mean optical 
depth of our models was scaled to agree with observations, but this normalization is in some dis- 
pute. However, changing this normalization has little effect on the median of the 6-distribution. For 
example, using sCDM (ag = 0.7), an increase in Tq from 0.225 to 0.35 at z = 3 causes the median 
b value to decrease from 20.8 km/s to 20.1 km/s, a change of only ~ 1 km/s. One possibility might 
be to change the ionization history of the universe such that pressure broadening would widen the 
absorbing structures. Another might be to change the density structure through the power spec- 
trum of the cosmology itself. However, with the suite of models considered here it seems difficult 
for a single model to give good agreement with both the column density and 6-parameter data. 



6. Flux Statistics for Helium II 

Previous work (Zhang et al. 1997,1998; Croft et al. 1997b) indicates that He ii Lya (304 A) ab- 
sorption may be significant in regions where H i Lya is not. Thus the study of He ii Lya absorption 
in quasar spectra provides a unique probe of structure in the lowest density regions of the universe. 
Comparison of both H i and He ii absorption within the context of a given cosmological model may 
also yield important information about the spectral shape of the metagalactic UV radiation field 
and its redshift evolution. While current observations still struggle to obtain sufficient resolution 
to detect any but the broadest individual He ii Lya lines, it is still possible to determine mean 
statistics of the He ii flux and optical depth which are not so sensitive to instrumental resolution. 
We define the mean optical depth thsII = — ln(F) where {F) is the mean transmitted fiux (F = 1 
signifying complete transmission). In Figure ^ we present th^jj as a function of redshift for the 
sCDM models with varying power normalizations (top) and for the models in the model comparison 
study (bottom). Several trends are apparent. First all models produce a rapid rise in mean optical 
depth with increasing redshift (roughly a factor two between z = 2 and z = 3), with tCDM and 
CHDM rising slightly more steeply. This is consistent with previous work on a smaller number 
of hierarchical cosmologies (Zhang et al. 1997; Croft et al. 1997b) and with the interpretation 
that the observed optical depth is due primarily to absorption by gas in underdense regions. The 
redshift evolution of the optical depth is thus dominated by the change in the gas density due to 
universal expansion and (to a lesser degree) by the shape of the UV metagalactic ionizing back- 
ground (here assumed to be that of Haardt-Madau (1996) with frequency dependence oc z^~^'®). 
Second, for a given redshift z, models with less power on small fiuctuation scales have progressively 
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larger optical depths. This is again consistent with the interpretation that the absorption is due to 
gas in predominantly underdense regions since less gas in these low power models will have turned 
around and collapsed. 

The first observation of a flux decrement at the wavelength where the He ii Lya absorption 
should occur was made by Jakobson et al. (1994) using the HST Faint Object Camera to ob- 
serve quasar Q0302-003. They obtained a 90% confidence lower limit of fHeii > 1-7 at z = 3.286. 
Subsequently improved measurements using spectra from this same quasar were made by Hogan 
et al. (1996) with the Goddard High Resolution Spectrograph on the HST and by Heap et al. (1998) 
using the Space Telescope Imaging Specrometer (STIS). STIS provides better sensitivity and 
background determinations than previous measurements. Davidsen, Kriss & Zheng (1996) used 
the Hopkins Ultraviolet Telescope to study the average He ii opacity in the spectrum of quasar 
HS1700-I-64 over the redshift interval 2.2 < z < 2.6 (lower than that available with HST). They 
find THeii = 1-0 ± 0.07, although as shown in Figure 3 of Croft et al. (1997b) there is considerable 
scatter when the wavelength range is divided into lOA bins. Measurements of He ii absorption have 
also been made by Anderson et al. (1998) with STIS using the spectrum of quasar PKS 1935-692. 
Although the number of lines of sight studied so far are limited and thus a detailed comparison of 
observations with our model simulations (that average over 300 lines of sight) is premature, these 
data are also presented in Figure ^. On face value these data favor higher optical depths and thus 
models with lower fluctuation power. However, none of the simulation models presented here can 
reproduce the apparent break at z = 3 in the optical depth observed by Heap et al. (1998). If this 
break persists it would most likely signal a departure from the Haardt-Madau quasar reionization 
spectrum assumed here. 

For completeness and comparison with previous work (Croft et al. 1997b) we show in Figure 
^ (from top to bottom) the He ii Lya flux probability distribution functions at z = 4, 3, 2, respec- 
tively, for the models in the model comparison study. The distribution functions are calculated from 
the flux smoothed with window functions corresponding to the same Full Width at Half Maximum 
(FWHM) as STIS with high, 50 km s~^ resolution (left column) and low, 500 km s~^ resolution 



(right column). As Figure 21 shows, the shape of the flux PDF is highly dependent on the smooth- 
ing. We see however, that models with less fluctuation power on small scales have far fewer truly 
transparent regions (pixels with F near 1). For the high Av = 50 kms^^ resolution case, all models 
converge in the fully saturated regime {F < 0.05) as expected. 



7. Summary 

We have performed several simulations of the Lya forest using different background cosmologi- 
cal models, numerical codes and grid resolutions. Five different cosmological models were considered 
here: the standard flat critical density cold dark matter model (sCDM), a flat CDM model with 
a nonzero cosmological constant (ACDM ), an open CDM model (OCDM), a flat critical density 
CDM model with a tilted power spectrum matching both the COBE amplitude and small scale clus- 
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tering contraints (tCDM), and a flat critical density mixed dark matter model (CHDM). The high 
resolution shock capturing code Kronos was used with grid resolution Ax = 37.5 kpc (Az = 75 kpc 
for CHDM) for the benchmark calculations presented in this paper. Three of the models (sCDM, 
ACDM , OCDM with identical parameters) were also evolved with the artificial viscosity based 
code Hercules at the lower grid resolution Ax = 75 kpc. Both simulation techniques give similar 
results for statistics, such as the slope of the column density distribution, that are insensitive to 
grid resolution. 

We have presented results from several statistical analyses of absorption features present in the 
Lya spectra, both from the unprocessed optical depth data and from the reduced line lists. Explic- 
itly we have considered the optical depth and transmitted flux probability distribution functions, 
the cumulative optical depth distributions, the H i column density distributions, line number evo- 
lution, Doppler b parameter distribution, the average flux difference as a function of velocity (first 
moment of the 2-point flux distribution function), and the mean optical depth and flux probability 
distribution functions for He ii absorption. We find: 

1. Simple scaling laws describe the redshift evolution of the optical depth, flux PDF, the 
H I column density distribution and, in conjunction with the slope of the column density 
distribution, the line number evolution remarkably well. This demonstrates that most of the 
evolution of the Lya forest is a direct consequence of universal expansion. 

2. The shape of the optical depth PDF is strongly correlated to the amplitude of the density 
fluctuation spectrum. Differences between models my be significant in the observationally 
accessible region 0.02 < r < 4. Thus this statistic may be a useful discriminator among 
models. Similar conclusions hold for the related flux PDF. 

3. Cumulative opacity distributions for the models are strongly clustered at low optical depth 
thresholds and high z. Significant differences do occur for optical depth thresholds tq > 1, 
but these may be more difficult to observe. 

4. The column density distribution function is a robust statistic relatively insensitive to grid 
resolution and numerical method. Its redshift evolution is described well by the same naive 
scaling law that describes the evolution of the optical depth. The slope of the column den- 
sity distribution is sensitive to the amplitude of the power spectrum on scales roughly the 
size of the absorbers (~ 100 kpc). Models with less power at these scales produce steeper 
distributions in qualitative agreement with semi- analytic arguments (Hui, Gnedin & Zhang 
1997). A preliminary comparison with data favors models with more power (sCDM, ACDM 
, OCDM) over those with less power (tCDM, CHDM). 

5. All models show comparable evolution for the number of lines above a given H I column 
density threshold in reasonable agreement with the data. Thus this statistic is not a sensitive 
discriminator among models. 
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6. Although the shape of the Doppler b parameter distribution is well reproduced by all the 
models, the median of the distribution for sCDM, ACDM , and OCDM models is well below 
observed values. The median of the b parameter increases for models with less power on small 
scales. Thus the observations favor low power models, such as tCDM, making it difficult for 
any model considered in this study to simultaneously give good agreement with both the 
H I column density and b parameter data. This discrepancy is confirmed as well in the 
nonparametric first moment of the two-point flux distribution function and so is not solely 
the result of the line- fitting algorithm employed. The solution to this problem may require 
modification of the reionization history of the universe to produce more pressure broadening of 
the absorbing structures or a modification of the power spectrum of the underlying cosmology 
itself. 

7. All of the models simulated in this study produce a rapid rise in He ii mean optical depth with 
increasing redshift consistent with the interpretation by previous work (Zhang et al. 1997, 
Croft et al. 1997b) that the observed optical depth is due to absorption by gas in underdense 
regions where universal expansion dominates the evolution of the gas density. Models with less 
power on small scales (tCDM, CHDM) produce larger mean He ii optical depths. Preliminary 
comparison with the data tends to favor these low power models. However, none of the models 
can reproduce the break seen by Heap et al. near z = 3. If this break persists in the data, it 
would most likely reflect that the Haardt- Madau (1996) form for the metagalactic UV ionizing 
background, based on homogeneous reionization by quasars alone in a clumpy medium, must 
be modified. 
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Fig. 1. — Power spectra for the models in the model comparison study (sCDM ug = 0.7, ACDM , 
OCDM, tCDM, and CHDM) plotted for a hnear evolution to z = 3. 
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Fig. 2. — Effective HI optical depth as a function of redshift for simulations of the model comparison 
study normalized to f = 0.3 at z = 3 compared to recent data by Kirkman & Tytler 1997 (denoted 
by stars). 
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Fig. 3. — (top) Optical depth probability distribution function for the high spatial resolution sCDM 
model with ug = 0.7 for redshifts z = 2,3 and 4. (middle) The same optical depth probability dis- 
tribution functions with r scaled by f{3)/f{z) where f{z) = (1 + z)'^'^ / (T (z) H) is the naive scaling 
predicted assuming no redshift evolution of the gas overdensity. (bottom) The same distributions 
with additional scaling due to the z evolution of the peak of the density distribution shown in 
Figure |^. 
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Fig. 4. — Redshift evolution of the gas overdensity for the high spatial resolution sCDM model 
with as = 0.7. 
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Fig. 5. — Predictions for rdP/dr for those cosmologies in the model comparison study (top). 
Dependence of rdP/dr on the fluctuation power spectrum in sCDM models with varying ug and 
simulation grid resolution (bottom). 



10 









o 
/ 






q/ 


/ 




/o 


/o 


- 


/ 






- 











Fig. 6. — Correlation of 0",-, a measure of the width of the rdP/dr distribution, with (T34, the 
amplitude of the linear power spectrum on small scales as given in Equation [l|. 
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Flux 



Fig. 7. — (top)Flux probability distribution for the high spatial resolution sCDM model with 
o"8 = 0.7 for z = 2,3 and 4. (bottom) The same sCDM distributions as above where both the 
distributions and the flux are scaled to the z = 3 values predicted by the simple scaling of r given 



by Equation 12 
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Fig. 8. — (top)Flux probability distributions at z = 3 for sCDM with spatial resolution Ax = 
37.5 kpc and ag = 0.7 and 0.3, respectively, (bottom) Flux probability distributions at z = 3 for 
the models in the model comparison study. 
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Fig. 9. — Scatter plot of the line center opacity as a function of column density for the sCDM 
model run with the Hercules code. Also shown (solid line) is the average line center opacity in the 
24 density bins of width Alog Nhi = 0.24 cm~^ . The cutoff in opacity at Tc > 0.05 is due to the 
spectral threshold e^"^ = 0.95 used in identifying the absorption features. 
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Fig. 10. — Fraction of pixels greater than three threshold optical depths, tq = 0.1 (circles), tq = 1.0 
(squares), and tq = 7.0 (triangles), plotted against redshift. The curves for tCDM and CHDM are 
indistiguishable for the lowest two thresholds. 
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Fig. 11. — Fraction of pixels exceeding a threshold optical depth tq at z = 3 for models in the 
model comparison study. The curves for tCDM and CHDM are indistinguishable. 
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Fig. 12. — (top) Uncut H l column density distributions for the high resolution sCDM erg = 0.7 
model for redshifts z = 2, 3 and 4; (bottom) The same H i column density distribtions with A'^hi 



scaled according to Equation |15| to z = 3 
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Fig. 13. — (top) H I column density distribution for high spatial resolution (Ax = 37.5 kpc) sCDM 
models with o^ = 0.3 (dashed) and ag = 0.7 (solid), (bottom) H i column density distributions at 
redshift z = 3 for (Kronos) simulated spectra in the model comparison study. Observational data 
are from Kirkman & Tytler (1997) and Kim et al. (1997). 
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Fig. 14. — Evolution of the number of lines with column densities greater than log A^hi > 13(top), 
logiYm > 13.5 (middle), logiVm > 14 (bottom) per unit redshift for the sCDM o-g = 0.7 (solid), 
ACDM (dashed), OCDM (dotted), tCDM (dot-dashed), and CHDM (long dashed) models of the 
model comparison study. The numerical results are compared against the observed data of Kulkarni 
et al. 1996 (filled circle), Hu et al. 1995 (filled square) and Lu et al. 1997 (filled triangle). 
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Fig. 15. — (top) Scaling predictions (lines) for the line number evolution with redshift for hnes 
with column densities greater than logA'^Hi > 13, logA^Hi > 13.5, logiVni > 14 per unit redshift 
computed for the high spatial resolution sCDM a% = 0.7 model and compared to the simulation 
data (symbols). The scaling prediction is normalized to the simulation data at z = 3. (bottom) 
Scaling predictions (lines) for the line number evolution with redshift for lines with column densities 
greater than logiVni > 14 for the models in the model comparison study compared to simulation 
data (symbols) for the same models. 
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Fig. 16. — Scaling law predictions of the redshift evolution of the metagalactic UV radiation field 
for models in the model comparison study compared to that of the Haardt-Madau field used in 
the simulations. The amplitude of the predicted radiation field is normalized to agree with the 
Haardt-Madau field at z = 3. 
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Fig. 17. — Doppler parameter distributions at redshift z = 3 from the high resolution (Ax = 
37.5 kpc) simulated data, (top) Doppler parameter distributions for sCDM, ACDM , OCDM, 
tCDM models of the model comparison study. The histogram is data from Kim et al. (1997) for 
z = 3.35. (bottom) Doppler parameter distributions for sCDM us = 0.3 and 0.7 to investigate the 
dependence on the fluctuation power spectrum. 
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Fig. 18. — Evolution of the median Doppler parameters for lines with column densities 10 

A^Hi < 10^^ (top) and 10^^'^ < A^hi < 10^"^ (bottom) for the high resolution models of the model 

comparison study. Thick solid lines represent data from Kim et al. (1997). 
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Fig. 19. — The average flux difference as a function of velocity for pixels in the flux range (from top 
to bottom in each panel) 0.0 < F < 0.1, 0.3 < F < 0.6, and 0.6 < F < 1.0 for the high resolution 
sCDM cjg = 0.7 and erg = 0.3 models (upper panel) and the high resolution models of the model 
comparison study (bottom panel). 



10 



10" 



10"' 



10" 



10"' 



Ax=75 kpc, Oj=0.7 

Ax=75 kpc, Oj=0.5 

Ax=75 kpc, Oj=0.3 

♦ HAR 
■ DKZ 
A Anderson 

* Heap 



— sCDM, a^= 

- ACDM 
OCDM 

- tCDM 

- CHDM 

♦ HAR 
■ DKZ 

▲ Anderson 

* Heap 





1.0 



2.0 



3.0 



4.0 



Fig. 20. — Evolution of the mean optical depth fueii for He ii Lya absorption for the sCDM 
models with cg = 0.7, 0.5, and 0.3, respectively (top) and for the models of the model comparison 
study (bottom). Also shown are observational data with error bars from Hogan et al. (HAR, 
fllled diamond), Davidson et al. (DKZ, fllled square), Anderson et al. (fllled triangle), and Heap 
et al. (stars). 
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Fig. 21. — He II Lya flux probability distribution function smoothed with window functions of 50 
kms"^ (left) and 500 knis~^ (right), respectively for redshifts (from top to bottom) z = 4, 3, and 2 
for models in the model comparison study. 
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Qo 


Qa 


qo 
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0"8/j-i 


Ax (kpc) 


17b /i2 


0"34 


sCDM 


1 





0.5 


0.06 


0.5 


1 


0.7 


37.5 


0.015 


1.89 


ACDM 


0.4 


0.6 


-0.4 


0.0355 


0.65 


1 


1.0 


37.5 


0.015 


2.03 


OCDM 


0.4 





0.2 


0.0355 


0.65 


1 


1.0 


37.5 


0.015 


2.50 


tCDM 


1 





0.5 


0.07 


0.6 


0.81 


0.5 


37.5 


0.025 


1.09 


CHDM 


1 





0.5 


0.07 


0.6 


0.98 


0.7 


75 


0.025 


1.14 



Table 1: Physical parameters of the different cosmological models. 0,q is the total density parameter, 
r^A = A/3Hq the cosmological constant density parameter, qq = Q.q/2 — JIa the deceleration 
parameter, 0.^ the baryonic mass fraction, h the Hubble parameter, n the slope of the primordial 
density perturbation power spectrum, cxg/j-i the fluctuation normalization in a sphere of radius 8/i^^ 
Mpc, Ax the comoving simulation spatial resolution in units of kpc, Of,/i^ is the baryon density in 
physical units (independent of Hq) and the last column is proportional to Gnedin's (1998) measure 
of power at small scales. For CHDM the fraction of the energy density carried in neutrinos is 0.2. 
The sCDM, ACDM and OCDM models were also simulated with the Hercules codes at lower spatial 
resolution (Ax = 75 kpc). 



Model 


Z 


di 


/?HGZ 


/3h 


/3f 


/3i 


/32 


sCDM 


2 


1.64 ± 0.03 


1.49 


1.84 ±0.23 


1.87 ±0.05 


1.63 ± 0.04 


1.98 ± 0.08 




3 


1.56 ±0.03 


1.55 


1.77 ±0.10 


1.71 ±0.03 


1.52 ±0.04 


1.79 ± 0.05 


ACDM 


2 


1.61 ±0.04 


1.50 


1.80 ±0.17 


1.74 ±0.04 


1.56 ±0.05 


1.73 ± 0.08 




3 


1.48 ± 0.04 


1.56 


1.78 ±0.12 


1.66 ±0.04 


1.42 ±0.03 


1.74 ±0.03 


OCDM 


2 


1.56 ±0.04 


1.47 


1.78 ±0.17 


1.71 ±0.04 


1.50 ±0.03 


1.78 ± 0.08 




3 


1.41 ±0.03 


1.52 


1.45 ±0.10 


1.58 ± 0.04 


1.41 ±0.05 


1.71 ± 0.04 


tCDM 


2 


1.89 ± 0.04 


1.63 


1.93 ±0.18 


1.90 ± 0.05 


1.86 ± 0.06 


1.83 ±0.13 




3 


1.72 ± 0.08 


1.70 


2.03 ±0.12 


1.95 ±0.06 


1.64 ±0.13 


2.03 ±0.10 


CHDM 


2 


1.87 ±0.04 


1.66 


1.88 ±0.17 


1.95 ±0.05 


1.85 ±0.05 


1.90 ±0.14 




3 


1.65 ± 0.08 


1.67 


1.89 ±0.10 


1.92 ±0.06 


1.57 ±0.13 


2.01 ± 0.08 


Kim et al. 


2.31 


1.35 ± 0.03 














2.85 


1.39 ±0.26 




1.46 ± 0.07 


1.46 








3.35 


1.59 ±0.13 














3.7 








1.55 







Table 2: Determinations of the H i column density distribution slope over various column density 
ranges: (3i over the range 10^^-^ < A^'hi < 10^^-^ cm~^ , /?h over the range 10^^-^ < A'^hi < 10^"'-^ 
cm~^ , /3f over the range 10^^'^ < A^hi < 10"*^^ cm^^ , /3i and (32 represent the slopes found by 



^1^ cm^^and 10^^ < A^m < 10 



16 



splitting the distributions into two halves, 10 < A'^hi < 10 
cm~^ , respectively. The column labeled /3hgz is the prediction of the semi-analytic model of Hui 
et al. (1997). The last row provides the measured values reported by Kim et al. (1997). All quoted 
errors are 2a. 
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Model 


713 


713.5 


714 


7h 


li 


sCDM 


2.48 ±0.13 


3.49 ±1.10 


4.27 ±1.65 


2.97 


2.29 


ACDM 


2.42 ±0.36 


3.38 ± 0.92 


4.09 ±1.42 


2.91 


2.02 


OCDM 


2.18 ±0.89 


2.96 ±0.78 


3.71 ±0.95 


2.73 


1.72 


tCDM 


2.60 ±0.79 


4.29 ±1.10 


5.30 ± 2.79 


3.49 


3.28 


CHDM 


2.32 ±0.97 


4.09 ± 0.62 


5.36 ±2.25 


3.79 


3.10 


combined 


1.61 ±0.66 


2.82 ±0.83 


3.47 ±1.23 






Kim et al. 








2.78 ±1.42 


1.29 ±0.29 



Table 3: Slope of the predicted line number evolution for models in the model comparison study as 
a function of redshift fit to the form dN/dz oc (1 ± zy for various fixed column density thresholds: 
iVm > 10^^ (713), ^Hi > 10^^-^ (713.5), Nm > 10^"' (714), for 2 < z < 4 and column density ranges: 
]^Qi3.77 ^ jsi^^ < ]^Qi6 (^^)^ ]^oi3-i < Nui < 10^'' {-/e) for z = 2 and 3. The same fits are presented for 
the combined line lists of Kulkarni et al. , Hu et al. , Kirkman &: Tytler, and Lu et al. (combined) 
and for data from Kim et al. (1997). All quoted errors are 2a. 



